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Abstract 

Context. The interaction between stellar winds and the interstellar medium (ISM) can create complex bow shocks. The photometers 
on board the Herschel Space Observatory are ideally suited to studying the morphologies of these bow shocks. 
Aims. We aim to study the circumstellar environment and wind-ISM interaction of the nearest red supergiant, Betelgeuse. 
Methods. Herschel PACS images at 70, 100, and 160/im and SPIRE images at 250, 350, and 500 /im were obtained by scanning 
the region around Betelgeuse. These data were complemented with ultraviolet GALEX data, near-infrared WISE data, and radio 
21 cm GALFA-HI data. The observational properties of the bow shock structure were deduced from the data and compared with 
hydrodynamical simulations. 

Results. The infrared Herschel images of the environment around Betelgeuse are spectacular, showing the occurrence of multiple arcs 
at ~6-7' from the central target and the presence of a linear bar at ~9'. Remarkably, no large-scale instabilities are seen in the outer 
arcs and linear bar. The dust temperature in the outer arcs varies between 40 and 140 K, with the linear bar having the same colour 
temperature as the arcs. The inner envelope shows clear evidence of a non-homogeneous clumpy structure (beyond 15"), probably 
related to the giant convection cells of the outer atmosphere. The non-homogeneous distribution of the material even persists until the 
collision with the ISM. A strong variation in brightness of the inner clumps at a radius of ~2' suggests a drastic change in mean gas 
and dust density ~32 000yr ago. Using hydrodynamical simulations, we try to explain the observed morphology of the bow shock 
around Betelgeuse. 

Conclusions. Different hypotheses, based on observational and theoretical constraints, are formulated to explain the origin of the 
multiple arcs and the linear bar and the fact that no large-scale instabilities are visible in the bow shock region. We infer that the 
two main ingredients for explaining these phenomena are a non-homogeneous mass-loss process and the influence of the Galactic 
magnetic field. The hydrodynamical simulations show that a warm interstellar medium, reflecting a warm neutral or partially ionized 
medium, or a higher temperature in the shocked wind also prevent the growth of strong instabilities. The linear bar is probably an 
interstellar structure illuminated by Betelgeuse itself. 

Key words. Stars: AGB and post-AGB, Stars: mass loss, Stars: circumstellar matter, Stars: individual: Betelgeuse 



1. Introduction 

For decades, it has been understood that the mass lost by red 
giants and supergiants dominates the mass return to the Galaxy 



( Maeder|1992| >, yet the actual mechanisms of such stellar mass 
loss have remained a mystery. The generally accepted idea on the 
mass-loss mechanism for asymptotic giant branch (AGB) stars 
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* Herschel is an ESA space observatory with science instruments pro- 
vided by European-led Principal Investigator consortia and with impor- 
tant participation from NASA. 



is based on pulsations and radiation pressure on newly formed 
dust grains. However, oxygen-rich AGB stars suffer from the 
so-called 'acceleration deficit' dilemma, which states that mass- 
loss rates due to the formation of silicate dust alone are orders 
of magnitude less than observed ones (Woitke 2006b). The for- 
mation of both carbon and silicate grains (Hofner & Andersen 
|2007|) or micron-sized Fe-free silicates (Hofner 2008; Norris 
|et al. 2012 1 have been proposed as possible solutions to this 
dilemma. The proposed models for the mass loss of AGB stars 
are very unlikely to be applicable to red supergiants (RSGs) 
since they are irregular variables with small amplitudes. 
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Processes linked to convection, chromospheric activity, or 
rotation might play an important role as triggers for the mass 
loss in RSGs. Josselin & Plez (2007 ) proposes that turbulent 
pressure generated by convective motions, combined with radia- 
tive pressure on molecular lines, might initiate mass loss, but 
Alfven waves generated by a magnetic field might also con- 
tribute (Hartmann & Avrett[T9 84). The last few years have pro- 
vided more and more evidence that, at least for a fraction of 
evolved stars, the mass is not lost in a homogeneous way, but 
irregular and clumpy structures prevail (e.g. Wei gelt et al.|2 002 ; 
Kervella et al.|[2009 Decin et "aT7||201 1| >. This might have seri- 



ous consequences for the total amount of mass lost during these 
late evolutionary stages. Solving the riddle of the mass loss of 
these (higher mass) targets is important in the framework of stel- 
lar evolution as a key to estimating when a supernova explosion 
may occur. 

decade (as CRIRES, IRAM, APEX etc.) made it possible to 
study these chemical processes in detail, and to refine our es- 
timates of the total chemical enrichment of the ISM by these 
evolved stars. 

When these evolved stars with their surrounding winds 
move through the interstellar medium (ISM), a wind-ISM in- 
terface structure is expected in the form of a bow shock, a 
kind of cometary like structure pointing in the direction of mo- 
tion. Photometric observations performed with the IRAS tele- 
scope (Noriega "-Crespo et al.|1997) and la ter confirmed with the 
Spitzer Space Telescope ( |Ueta et al . 2006 2008 ) and the Galaxy 
Evolution Explorer satellite (GALEX Marti n et al.|2007| [Sahai 
|& Chrorlo poulos 2010) indeed showed the presence of bow 
shock structures around evolved stars. However, these early ob- 
servations of the far-IR and UV-bright emission structures lacked 
any spectral resolution, and since they only had poor spatial res- 
olution, allowed only basic morphological studies. 

Thanks to the capabilities of the instruments on board the 
Herschel Space Observatory (Pil bratt et al.||2010) , one is now 
able to image the circumstellar envelopes (CSEs) and the CSM- 
ISM interaction regions with unprecedented detail. In the MESS 
GTKP ( |Groene wegen et al. 201 1| > several known and new in- 
teraction objects have been imaged with the PACS and SPIRE 
instruments (e.g . CWLeo, |Ladjal et al.|2010)|Decin et al.|2011| 



Cox et aL|2012i, among which several supergiants. Betelgeuse 



(a Ori, M2 lab) is the closest supergiant, and it serves as the 
prototype of the red oxygen-rich supergiants. It has been an im- 
portant astrophysical laboratory for many decades and its enve- 
lope has been studied with a variety of observational techniques. 
Betelgeuse is, therefore, a prime target for Herschel to study the 
still unknown mechanisms that lead to mass loss and eventually 
result in an energetic bow shock when the wind collides with the 
ISM. 

1.1. Betelgeuse 

Many detailed studies have already described some characteris- 
tics of the complex atmosphere, chromosphere, and dusty enve- 
lope of Betelgeuse. In this section, a summary is given of those 
properties that are important for this study. 

Betelgeuse is a low-amplitude variable, whose variation was 
first described in 184 by Sir John Herschel ( jHerschellfl 840| l . 
Schwarzschild ( 1975| ) attributed these fluctuations to a changing 
granulations pattern formed by a few giant convection cells cov- 
ering the surface of the star. Recently, a magnetic field of about 1 
Gauss has been detected using the Zeeman effect ( |Auriere et al. 
2010). It is very faint, but may play an important role in shaping 
the convection of Betelgeuse. 



Using NRAO's Very Large Array (VLA), combined with 
Hipparcos data, Harper et al. ( 2008| ) determined the parallax of 
Betelgeuse to be 5.07±1.10mas (or a distance of 197±45pc). 
The most likely star-formation scenario for Betelgeuse is that it 
is a runaway star from the Orion OBI assocation and was origi- 
nally a member of a high-mass multiple system within Ori OB 1 a 
( |Harper et al.|2 008 ). The galactic space motion in the local stan- 
dard of rest (LSR) is U' = -12.7, V' = +0.5, W = +28.2 km/s. 
Betelgeuse was probably formed about 10 to 12 million years 
ago from the molecular clouds observed in Orion, but has 
evolved rapidly owing to its high mass. No consensus has yet 
emerged regarding the star's mass, but most studies report val- 
ues betwee n 1OM and 20 M Q (e.g., |Dolan et al.|2008"}|Neilson| 
|etal|2011D . 

Ultraviolet and Ha images prove there is a chromosphere 
with a spatial extent of ~4 R^ ( |Altenhoff et al.|1979||GoId berg 
|1981 Hebden et al.|1987|), where the hot plas ma has a temper- 



ature of around 5000 K (jGilliland & Dupree 



1996). However, 

observations made with the VLA by |Lim et al. (19981) suggest 
that much cooler (^1000-3000 K) gas co-exists at the same ra- 
dial distances from the central target, and must be much more 
abundant because it dominates the radio emission. In addition, 
narrow-slit spectroscopy of Betelgeuse at 11.15 /xm with the ISI 
by |Verhoelst et al.] ( |2006[ ) reveals that silicate dust forms at dis- 
tances beyond ~20R*, but that AI2O3 may form as close as 
^2R*. This apparent dichotomy is solved when adopting a sce- 
nario in which a few inhomogeneously distributed large con- 
vective cells are responsible for lifting the cooler photospheric 
gas into the atmospheres and for producing shock waves, which 
could heat localized, but spatially distributed regions of the at- 
mosphere to chromospheric temperatures ( |Lim et al.|1998 1. 

Observations taken over the past few years have proven that 
the regions just above the stellar surface show some clear devi- 
ations from spherical symmetry. IKervella et al. ( 2009 ) obtained 



adaptive optics images of Betelgeuse in the 1.04-2.17 /im range 
in ten narrow-band filters. These images show that the CSE of 
Betelgeuse has a complex and irregular structure, with in par- 
ticular a bright 'plume' extending in the southwestern quadrant 
up to a radius of at least six times the photosphere. They pro- 
posed that the 'plume' might be linked to either the presence 
of a convective hot spot on the photosphere or to polar mass 
loss possibly due to rotation of the star. Using AMBER on the 
VLTI, |Ohnaka et aT] d2009) ) found from CO lines that the at- 
mosphere has an inhomogeneous velocity field, which could by 
explained by the gas moving outward and inward in a patchy pat- 
tern. This gives support to the idea of convective motion in the 
upper atmosphere or intermittent mass ejections in clumps or 
arcs. Recently, Ker vella et al.| ( [20TT] ) observed Betelgeuse with 
VLT/VISIR in the N and Q band. The images showed a bright, 
extended and complex CSE within ^2.5" of the central star with 
the appearance of a partial circular shell between 0.5 and 1.0", 
several knots, and filamentary structures. The ring-like structure 
might be related to the dust condensation zone. 

Although some dust species, such as silicates and corundum, 
have been identified in the CSE around Betelgeuse, the envelope 



shows a remarkable deficiency in dust ( |Skinner & Whitmore 
|1987[ [Verhoelst et al.| [2006). In addition, the molecular con- 
tent is unusually low. There is evidence that only a fraction of 
carbon is associated into CO (Huggins |1987[ ). The extended, 
warm chromosphere is clearly unfavourable for the survival of 
molecules and the formation of dust species. Detailed studies by 



Glassgold & Huggins|(|T986) show that the most common form 



of all species are neutral atoms and first ions. 
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Figure 1: Herschel PACS and SPIRE images of Betelgeuse. North is up, east to the left. The field-of-view (FOV) is ~49'x29'. 



A parsec-size shell around Betelgeuse was detected with 
IR AS by|Noriega-Crespo et al.| ( fl997l a nd later on with AKARI 
by |Ueta et aT ( 2008 1 (see Appendix [A}. The shell is asymmet- 
ric with a mean radius of 6' and the outermost structure at ~7'. 



Noriega-Crespo et al. (19971) proposed that the shell is confined 
by the ram pressure of the ISM, hence represents the bow shock 
around Betelgeuse. To the north-east of Betelgeuse, a ridge or fil- 
ament was detected, which was said to belong to the surrounding 
ISM. Neither IRAS nor AKARI could clearly resolve the bow 
shock and filament structure. The derived dust temperatures in 
the bow shock and filament vary between ~10 K and 40 K ( |Ueta| 
et al.|2008) . 

Recently, |Le Bertre et al.| ( |2012| l has detected a detached H I 
gas shell of ~2' in radius surrounding Betelgeuse, and reported 
the detection of atomic hydrogen associated with the far-infrared 
arc located 6' north-east of Betelgeuse. 

In this paper, we focus on the morphological appearance of 
the inner envelope and bow shock surrounding Betelgeuse. The 
thermodynamic al structure and chemical content of the envelope 
and bow shock will be discussed in |Decin et al.| ( |2012| l (from 
here on called Paper II). In Sect. [2] the different observations 
with their respective calibrations are described. The observa- 
tional properties of the inner envelope and bow shock structure 
are analysed in Sect. [5] With the aid of hydrodynamical simula- 
tions, the origin of the multiple arcs is discussed in Sect. [4] The 
origin of the linear bar is described in Sect. [5] The conclusions 
are summarized in Sect. [6] 



2. Observations 

In this section, an overview is given of the different observations 
toward Betelgeuse and its extended envelope, the data-reduction 
for each data-set is described, and some first results are deduced 
from the data. 



2. 1. Herschel PACS and SPIRE images 

Infrared images were obtained using the Photodetector Array 
Camera and Spectrometer (PACS, Poglitsch 'et al.|2010] | and the 



Spectral and Photometric Imaging Receiv er (|Griffin et al.| 2010 
SPIRE) on board the Herschel satellite ( |Pilbratt et al.pOlO) , 
and they are part of the MESS guaranteed-time key programme 
( Groenewegen et al.|201 1) . For both instruments, the 'scan map' 
observing mode was used with medium (20"/s) scan speed in 
the PACS 70, 100, and 160 /xm filter settings and 30'7s for the 
SPIRE 250, 350, and 500 fim filters. To create a uniform cover- 
age, avoid striping artefacts, and increase redundancy, two ob- 
servations at orthogonal scan directions were concatenated. The 
PACS images were obtained on September 13, 2010 (OBSIDS 
1342204435 and 1342204436) and on March 29, 2012 (OBSIDS 
1342242656 and 1342242657), the SPIRE images on March 11, 
2010 (OBSID 1342192099). 

The data reduction was performed using the Herschel 
Interactive Processing Environment (HIPE) (see |Groenewegen| 



et al.| 201 1 for more details). To remove low-frequency noise 
which causes additive brightness drifts, the SCANAMORPHOS 
routine (Roussel 2012]) was applied. Deconvolution of the data 
was done using the method as described in Ottensamer et al. 



d20TT) . 

The PACS instrument offers a pixel size of 3.2" in the 70 
and 100 /xm bands and of 6.4"in the 160 zxm band, but the im- 
ages presented in this paper are oversampled by a factor of 3.2, 
resulting in a sampling of 1" and 2" per pixel, respectively. The 
PACS FWHM at 70, 100, and 160/xm is 5.7", 6.7", and 11.4", 
respectively. The pixel size of the three SPIRE bands is 6" at 
250 fim, 10" at 350 /xm, and 14" at 500 /xm, with an FWHM of 
the beam size being 18.1", 25.2", and 36.6", respectively. 

The individual scan maps are plotted in Fig. [T] and a false- 
colour image displaying all the PACS data is shown in Fig. [2] 
Remarkably, a multiple arc-like structure is detected at ^6-r 
away from the central star, together with an enigmatic linear bar- 
like structure at ~9'. A closer zoom into the arc-like structure is 
displayed in Fig. [3] The inner radius of the closest arc is at 280" 
(as measured from the central target). The other most prominent 
arcs have an inner radius at 280", 310", 350" and 375". The 
brightest parts of the arcs are reasonably well fitted by concentric 
ellipses (see Fig. |B.l| in the Appendix [B|. The width of the arcs 
is typically ~20 . Clearly, substructures of the size of several 
arcseconds, probably related to density differences, are seen in 
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Figure 2: Composite colour image of the Herschel PACS images of Betelgeuse. North is up, east to the left. Blue is the PACS 70 fim 
image, green is PACS 100 /im, and red PACS 160 /im. The black arrow indicates the direction of the space velocity of the star. The 
multiple arc-like structure is situated at ~6-7' to the north-east of the central target, the linear bar at ~9'. The contrast in the figure 
is best visible on screen. 



the arcs. The arc -like structure is probably related to the CSM- 
ISM interaction phase, when the circumstellar material collides 
at a wind speed of ~15 km/s with the ISM, creating a bow shock 
(see Sect.[4]for more details). The inner ridge of the linear bar is 
at 535" from the central star. The length is ^1600". The width is 
^30" in the northern part; moving down towards the south-east, 
the bar widens and splits into two parts. The angle between the 
direction of space velocity (drawn in Fig. |2|i and the linear bar is 
-80°. 



In the whole sample of 78 AGB stars and red supergiants sur- 
veyed in the framework of the MESS GTKP, Betelgeuse is the 
most spectacular source, because it is the only one showing this 
multiple-arc structure and the presence of a linear bar. The linear 



bar and the full arc were already detected with IRAS ( Noriega- 
|Crespo et al.|1997) and AKARI flUeta et al.|2008) , but their in- 
dividual structures could not be resolved (see Appendix |A|. 

To emphasize the shell morphology in the inner envelope, we 
removed the extended envelope halo by subtracting a smooth, 
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Figure 3: Zoom into the multiple arc-like structure seen in the 
bow shock region of Betelgeuse as detected with the PACS 
70 fim filter. North is up, east to the left. The FOV is 855" x 765". 
Flux units are in Jy/pixel. 



azimuthally aver aged profile repre sented by a power law r a 
for r <240" (see|Decin et al. 2.01 1}>. The method has been tested 



in detail to determine if artefacts from the point-spread function 
(PSF) still affect the data (Royer et al., in prep.). Only the region 
within 15-20" from the central target should be interpreted with 
care owing to the complex PSF and the bright contribution from 
the central target at these infrared wavelengths. The resulting 
image is shown in Fig. [4] The inner envelope structure is clearly 
non-homogeneous and testifies to a turbulent mass-loss history. 
A detailed discussion of the inner envelope structure is presented 
in Sect. O 

2.2. WISE data 

The complex environment around Betelgeuse has also been im- 
aged by the Wide-field Infrared Survey Explorer (WISE, Wright 
et al.||2010 l. WISE has four photometric bands centred at 3.4, 
4.6, 12, and 22 /an, with an angular resolution of 6.1", 6.4", 
6.5", and 12", respectively. The 5er point source sensitivity in 
unconfused regions on the ecliptic is better than 0.08, 0.11, 
1, and 6mJy, respectively. The WISE All-Sky Data were re- 
leased on March 14, 2012. The WISE data images are shown 
in Fig. [5] The extreme brightness of Betelgeuse at these infrared 
wavelengths makes it very difficult to find any nearby associ- 
ated emission structures. The linear bar is visible only at 12 
and 22 /im, with a (background-subtracted) surface brightness 
of — 1.2MJy/sr and — 1.9 MJy/sr, respectively, and 5.2269e-5 
for the 22 ^m band. 



2.3. GALEX observations 

Betelgeuse was also observed with the Galaxy Evolution 
Explorer satellite (GALEX) on November 29, 2008. The 



pipeline-calibrated FUV and NUV images were retrieved from 
the GALEX archive. The bandpass is 1344-1786 A for an angu- 
lar resolution of 4.5" for the FUV filter, and 1771-2831 A for a 
6" angular resolution for the NUV filter (Morris sey et al.|2 005). 
The pixel size is 1.5" x 1.5". The integration time for Betelgeuse 
was 48916". Some very faint extra extinction accompanying the 
outermost arc around Betelgeuse was detected in the FUV filter 
(see Fig.|6]l. The FUV flux of the faint extra extinction is ^7% 
lower than in the adjacent pixels. The NUV filter is unfortunately 
satured by a ghost of the central star. 

Until now, bright UV emission in the wind-ISM interaction 



region has only been detected in two AGB stars: o Cet ( Martin 
|et al |2007| > and CW Leo ( |Sahai & Chronopoulos|20T0] l; a Ori 
is the first target for which extra UV extinction in the wind- 
ISM interaction zone is reported. No UV emission or extinc- 
tion is seen in the region of the linear bar. The origin of the 
FUV/NUV emission of the bow shock around o Cet has been 
attributed to H 2 molecules in the cold gas that are collision- 



ally excited by hot electrons from the post-shock gas (Martin 
|et al.|200"7] ). The origin of the FUV/NUV emission in the inter- 
action zone between the ISM and the stellar wind of CW Leo has 
not been modelled in detail, although it was suggested by Sahai 
|& Chronopoulos| ([2010) that also here collisionally excited H 2 
emission might be the origin. The FUV/NUV emission around 
CW Leo peaks at slightly larger distances from the central target 
than the PACS/SPIRE infrared flux excess (see Fig. 2 in Ladjal 
eTaT]|20T0] l. We note that for the planetary nebula NGC6720, 



dust and H 2 are co-spatial, and it is argued that H 2 has been 
formed on grain surfaces ( |van Hoof et al.|2010| l. If the FUV ex- 
tinction toward the bow shock around Betelgeuse is the FUV 
counterpart of the outer arcs as detected with PACS, this sug- 
gests that background UV radiation is absorbed by dust in the 
bow shock. 

2.4. GALFA Hi 21 cm observations 

Hi 21 cm observations in the vicinity of Betelgeuse were ex- 
tracted from the results of the Galactic Arecibo L-Band Feed 
Array H I (GALFA-H I) survey ( |Peek et al.|2"0TT) . GALFA-H I 
is a survey of the Galactic interstellar medium with a spatial res- 
olution of ~4', covering a large area (13 000 deg 2 ) with a high 
spectral resolution (0. 18km/s) from —700 < t>LSR < +700km/s 
in the 21 -cm line hyperfine transition of neutral hydrogen con- 
ducted at the Arecibo Observatory. ALFA is a seven-beam feed 
array, with the beams arranged on the sky in a hexagon. Each 
beam is slightly elliptical along the zenith-angle direction, ap- 
proximately 3.3 x 3.8' FWHM, with a gain of ^11 K/Jy for the 
central beam and ^8.5 K/Jy for the six other beams. Details on 



the data reduction are given in |Peek et aL 



2011). Typical noise 



cm/s channel, corre- 



levels are 80 mK RMS in an integrated 1 
sponding to an integration time of 30 sec. 

Several selected data cubes are displayed in Fig. [7] At the 
CO w LS r velocity of — 4km/s ( |De Beck et al.||2010| >, Hi emis- 
sion is clearly detected. In Fig. [8] we present the HI spec- 
trum obtained on the position of Betelgeuse. The spectrum is 
in good agreement with the Nangay RadioTelescope (NRT) and 
Leiden/Argentina/Bonn (LAB) data presented in Fig. 3 by |Le| 
|Bertre et al.| ( |2012| l. The emission peaks at 6.992 km/s, and 
no distinctive feature is seen in the velocity range of the cir- 
cumstellar CO emission. It illustrates that in the direction of 
Betelgeuse, the Hi emission is dominated by galactic emis- 
sion due to interstellar atomic hydrogen along the line of sight 
( |Le Bertre et al.||2012| ). The Hi spectrum was also inspected 
at different offset positions. Taking the central star in the 'on- 
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Figure 4: Oversampled PACS 70 /im image of Betelgeuse after subtraction of the smooth halo of the circumstellar envelope. North is 
up, east to the left. The pixel size is 1". In the left-hand panel, the power a has been taken as constant over the full region (a = 1 .78). 
The field-of-view (FOV) is 720" x 325". The red circle with a radius of 20" shows the region where some PSF artefacts are still 
visible; the green circle indicates a region with radius of 2'. In the right-hand panel two independent profiles were fitted before 
and beyond 50" to suppress the high flux values around 50" so as to better visualize the complex structure around 50-60" in the 
southeastern region. The field-of-view (FOV) is 970" x 640". Distances to a few asymmetric structures are indicated in red. 



PACS 70 mic 



WISE4.6MIC ' WISE 12MIC 





Figure 5: Comparison between the WISE observations and the PACS 70 /im image. An ellipse and line are added to guide the eye 
to the arc and bar-like structure. The pixel size of the WISE data is 1.375". The bright central star causes latents and glint artefacts 
in the images. 



position' and the 'off-positions' at different offsets gave the 
H I spectrum of Betelgeuse (see Fig. [9]). Although affected by 
interstellar confusion, emission at the CO UlSR is clearly de- 
tected in the spectra with the offsets taken in the east-west di- 
rection. Fig. [10] shows the H I spectrum at different offset po- 
sitions for which the off-position was taken in a region out- 
side the envelope, arcs or linear bar (i.e. at 11' to the west and 
10' to the north of the central target). The interstellar confu- 
sion is seen for instance in the emission feature at —3 km/s, and 
the absorption at +12 km/s. The Betelgeuse Hi line profile has 
an FWHM of ~3.5km/s. The integrated line intensity is ~2- 
7Jykm/s. Assuming a smooth outflow and using the standard 
relation, M Hl = 2.36 x 10~ 7 L> 2 / S v dv, with M H: in M Q , the 
distance D in pc, and J S v dv in Jykm/s, this translates to an 
envelope mass of ^0.02-0.07 M Q in atomic hydrogen at 197 pc. 
As discussed already, the density in the inner envelope is not 
homogeneous, for which reason this mass value should be con- 
sidered as very approximate. 



Significant emission is also seen at other velocities, of which 
a few examples are shown in Fig. [7] Using VLA data, |Le Bertre| 
|et al. (2012 1 note a seemingly spatial association between the 
H I emission integrated over the range —1.5 to 8.9 km/s and the 
IRAS 60 /im image (see their Fig. 9). This coincidence is not 
seen between the GALFA-H I and Herschel 70 /im data, proba- 
bly owing to the low spatial resolution of the GALFA-H I data 
(see upper panel in Fig.fTT). 

The bottom panels in Fig. [7] showing the emission at a ve- 
locity of 18 and 21. 17 km/s, notably display an alignment of the 
H I emission with the bar and arc-like structures detected in the 
Herschel PACS images (see also bottom panel in Fig. [TTJ. This 
H I emission probably has a galactic ISM origin, and might be 
part of cold, low-density, atomic gas structures (see also Sect. [5}. 
To gain further insight into the origin of the bar-like structure 
seen in the PACS observations, we compared the galactic Hi 
emission at different velocities in a larger region covering a 
right ascension from 05h39m56s to 05h59m33s and a declina- 
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Figure 6: Comparison between PACS 70 /im (left) and GALEX FUV image (right, after applying a boxcar smoothing with a kernel 
radius of 2). The white arrows indicate some faint extra extinction in the GALEX FUV image, which coincides with the outermost 
arc detected with Herschel. (The contrast in the figure is best viewed on screen.) 




Figure 7: Comparison between the GALFA-H I 21 cm observations (1' square pixels) and the 0.3 mJy countours of the PACS 70 fim 
image (red). The GALFA-H I data cube channel (CH) and corresponding velocity are indicated in the upper left corner. 



tion from 07°21'10.23" to 9°42'30.47". From these data, it is 
clear that the ISM contains a lot of small-scale structure and that 
most of the emission seen in the direction of Betelgeuse is not 
related in any way to the mass loss of Betelgeuse. Because it is 
situated beyond the Local Bubble, Betelgeuse probably lies near 
the edge of a cavity of low neutral hydrogen density and close 
to a steep gradient in ISM column density between 1= 180°- 
200° (see also Harper et al. 2008 1. The ISM small-scale structure 
in the direction of Betelgeuse has also been detected by Knapp 
|& Bowers] ( |1988| ) via the observation of discrete molecular CO 



clouds, probably associated with the A Ori molecular cloud com- 
plex. 



3. Observational properties of the infrared 
(extended) envelope around Betelgeuse 

In this section, some observational properties as deduced from 
the infrared data for the inner envelope structure, and the arcs 
and linear bar are described in Sects. [3~T|and|3.2| respectively. 
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Figure 9: a Ori spectra H I spectra obtained by taken the off-positions in different directions (as indicated in each panel). The offset 
position at 6' to the east is situated in the arc detected with Herschel. The right bottom panel shows the spectrum when the offset 
position is taken in the linear bar. The vertical line marks the CO i?lsr of Betelgeuse. 



3. 1. Inner envelope structure 

Figure [4] shows that the inner envelope structure of Betelgeuse is 
clearly non-homogeneous. The detection of non-spherical struc- 
tures has already been reported on a sub-arcsecond-scale (e.g. 
|Lim et al.|1998l|Kervella et al.|2009||20TT) (see Sect.[U). The 
Herschel images show the first evidence of a high degree of 
non-homogeneity of the material lost by the star beyond 15", 
which even persists until the material collides with the ISM. As 
described in Sect. |2.1| and shown in Fig. |4] very pronounced 
asymmetries are visible within 110" from the central star, al- 
though some weaker flux enhancements are visible until ^300" 
away. The typical angular extent for these arcs in the inner en- 
velope ranges from ^10°to 90°. Considering that asymmetries 
are seen on a subarcsecond-scale close to the stellar atmosphere 
and in the dust formation region, the Herschel images suggest 
that the arc-like structures seen in the free-flowing envelope of 
Betelgeuse (i.e. before the material collides with the ISM) are 
the relics of a non-homogeneous mass-loss process. 



The radial distance of the most pronounced asymmetries 
(within 2') agrees with the size of a ring of H I emission recently 
detected by |Le Bertre et aL ( 2012| l, although no counterpart to 
the H I emission plateau further away from the star in the south- 
western region is seen in the Herschel images. The Herschel and 
Hi observations seem to suggest that the mean gas and dust 
density (averaged over 360°) in the inner envelope drastically 
changed at ~2' (or 32000yr ago for an expansion velocity of 
3.5km/s). This might be due to, e.g., a change in mass-loss rate 
or the creation of an inner bow shock inside a fragmented, fila- 
mentary halo (see Sect. |4. 1 [ >. 

The angular extent of the arc-like structures in the inner en- 
velope of Betelgeuse is considerably smaller than for the AGB 
star CW Leo, the only other target in the MESS sample whose 



inner envelope has already been studied in detail (Decin et al. 
2011). In CW Leo, almost spherical, ring-like, non-concentric 
dust arcs where detected until 320", having an angular extent 
between ^40° and ^200°. The shells of CW Leo have a typical 
width of 5"-8", and the shell separation varies in the range of 
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Figure 10: H I spectra at different positions in the envelope, arcs, and linear bar. The spectrum of the 'off-position' taken at 11' to 
the west and 10' to the north of the central target has been subtracted. The position at 6' to the east is situated in the arc detected 
with Herschel. The right bottom panel shows the on-source Betelgeuse— off-source H I spectrum and the H I spectrum in the linear 
bar. The vertical line marks the CO «lsr of Betelgeuse. 
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Figure 8: GALFA-Hl spectrum obtained on the position of 
Betelgeuse. The vertical line marks the CO ulsr of Betelgeuse 
and the horizontal bar the total linewidth of the CO emission 
(twice the terminal velocity of the wind). 



formed non-isotropically, does not hold for Betelgeuse, which 
is an irregular variable with small amplitude variations. As de- 
scribed in Sect. 1 1.1 1 the mass-loss process might by convection- 
induced, which would yield locally strong variations in gas and 
dust density. 

Assuming that dust is the main contributor to the 
Herschel/PACS images (see Paper II), one can derive the dust 
temperature in the arcs from a modified blackbody of the form 
B U {T) ■ A"* 3 , as expected for a grain emissivity Q a b s ~ A~' 3 
with the emissivity index [3 ranging between oneQ (typical of, 
e.g., layered amorphous silicate grains, Kna pp et al.|[l 993) 



and two (typical for, e.g., crystalline silicate grains, 



Tielens & 



Allamandola| |1987| |Mennella et al.||1998) . When j3 is equal to 



2.0, the dust temperature for the clumps beyond 60 ranges from 
~25 to 65 K. For /3 equal to 1.0, the dust temperature increases 
to values between 35 K and 140 K (see Fig. 13 in Sect. 4.1 1. 



The Herschel images prove that clumpy structures are preva- 
lent over the full envelope and might eventually have an impact 
on the shape of the bow shock structure. 



~10"-35", corresponding to ~500-l 700 yr. This (again) sug- 
gests that the mass-loss mechanism of the AGB star CW Leo, 
based on pulsations and radiation pressure on dust grains that 



It is generally admitted from Kramers-Konig relations that 1 is a 



lower limit for the spectral index l Emerson 



19881. 
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Figure 1 1 : Comparison between the GALFA-HI emission, 
summed from -1.472 to 8.832 km/s (top) and from 17.848 to 
21.528 km/ (bottom), and the 0.3 mJy contours of the PACS 
70 /im image. 



3.2. Multiple arcs and linear bar 
3.2.1 . Orientation of the bow shock 



Using the Herschel PACS images, Cox et al. ( 2012| > derived a 
de-projected stand-off distance for the outer arcs of 4.98' and 
a position angle of 47.7°. Based on the distance of 131 pc and 
proper motion values of van Leeuwen ( 2007| l, a space motion 
vector inclination ilsr of 8.9° for a ^rad. LSR of 3.4 km/s was 
deduced. From the AKARI data, |Ueta etaT| ( |2008| l deduced a de- 
projected stand-off distance for the outer arcs of 4.8'±0. 1', a po- 
sition angle of 55°±2° and an inclination angle for the outer arcs 
of 56° ±4°, where the bow shock cone is oriented to the plane of 
the sky. From the difference between the galactic space-velocity 
components of the bow shock cone and the heliocentric galac- 



tic space-velocity components of Betelgeuse, |Ueta et al.| (2008 ) 
find that the ISM around Betelgeuse flows at ~1 1 km/s into the 
position angle of ~95° out of the plane of the sky ( towards us). 
Differences between the Herschel (Co x et al.|2012| l and AKARI 
( Ueta et al.|2008| l results are mainly because of using different 
values for the radial velocity, proper motion, and distance. Ueta 
et al.| ( [2008"l > and |Cox etaE] ( |20 1 2) have deduced that the pecu- 
liar velocity of the star with respect to the ISM is between 24 and 
33 km/s for an interstellar hydrogen nucleus density between 1.5 
and 1.9 cm~ 3 . 



3.2.2. Instabilities in the arcs 

The arc-like structure around Betelgeuse does not show clear 
evidence of large-scale instabilities, created by e.g. Kelvin- 
Helmholtz (KH), Rayleigh-Taylor (RT), Vishniac non-linear thin 
shell (NTSI Vishniac| [T994), or transverse acceleration (TAI 
Dga ni et al.||1996[ ) instabilities. This is in contrast to other tar- 
gets in the MESS sample as, e.g., R Leo for which the Herschel 
image shows clear signatures of RT instabilities that are slightly 
deformed by to the action of the KH instabilities (see Fig. |C.l 
in the Appendix |Cj. The upper limit on the maximum length 
for instabilities in the arcs around Betelgeuse as trace d from the 
Herschel images is ^30". Although suggested by Ueta et al. 
(2008), vortex shedding (Wareing et al. 2007b) downstream in 



the tail of the bow shock is not seen in the Herschel images. With 
a wind velocity of ^14.5 km/s (as deduced from low rotational 
CO transitions, see Paper II), the ratio between the peculiar ve- 
locity of the star, v+ and the wind velocity, v w , is more than 1 and 
instabilities in the shell of the bow shock are expected (Dgani 
et al.|1996| . This is also clear from the hydrodynamical simula- 
tions presented byjvan Marie et al.| ( [20ri) for parameters typical 
of a red supergiant resembling Betelgeuse. The non-occurrence 
of large-scale instabilities can put constraints on the physical 
parameters determining the morphology of the bow shock (see 
Sect.|42l. 



3.2.3. Dust temperature in the arcs and linear bar 

Both the multiple arcs and the linear bar are best visible in 
the PACS images. The surface brightness in the multiple-arc 
like structure and linear bar has been measured for the dif- 
ferent PACS filters using aperture photometry (see Table TJ. 
Uncertainties are ^15%, mainly dominated by the absolute flux 
uncertainty. The PACS and AKARI fluxes are consistent with 
each other for the arc-like region, with the exception of the 
AKARI 160 /im point which is a factor ^3 higher. However, 
the flux measurements in the bar-like structure are quite differ- 
ent (see Fig. 12 1. We note that both the arcs and the bar are ex- 
tremely faint in the AKARI LW images (see Appendix [Aj, and 
contamination cannot be excluded. These data points will there- 
fore be neglected. 

Neither [O I] line emission at 63 /im nor [C II] at 157 fjm is 
detected with Herschel in the bow shock region or linear bar (see 
Paper II for more details), which is the reason we assume that 
dust emission is the main contributor to the Herschel/PACS im- 
ages. Using a modified blackbody with a spectral index /3 vary- 
ing between one (typical for, e.g., layered amorphous silicate 



grains or amorphous carbon grains, Knap p et al.|| 1993 Dupac 
|et al.|2003] > and two (typical of, e.g., crystalline silicate grains or 
graphitic grains, Tielens & Allamand ola|| 1 987] [Mennella et al. 
|1998[ |Dupac et al.|2003 1, we derived the (mean) dust tempera- 
ture from the Herschel data using a Monte-Carlo uncertainty es- 
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Table 1 : Background-subtracted average surface brightness val- 
ues [in units of MJy/sr] as derived from the Herschel/PACS and 
WISE data (this paper) and from the AKARI as listed by Ueta 
et al. ( 2008| l for the arc and bar-like structure near Betelgeuse. 







ARC 


BAR 




70 fim 


17.9±2.69 


16.3±2.45 


PACS 


100 (im 


9.7±1.46 


9.8±1.47 




160 fim 


4.1±0.61 


5.3±0.80 




65 fim 


5-25 


~10 


AKARI 


90 (im 


5-20 


~10 




140 


5-10 


~15 




160 fim 


10-15 


~20 


WISE 


12 fim 




<0.74 


22 fim 




1.53±0.31 




Figure 12: Average surface brightness values derived from the 
Herschel/PACS, AKARI, and WISE data for the arc (left) and 
bar-like (right) structure around Betelgeuse. Shown in full black 
line is the best-fit modified blackbody to the Herschel data with 
a temperature of 85 K (/3= 1.0) for the arcs and 63 K (/3= 1.0) 
for the bar. The AKARI 140 and 160 fim data for the arcs and 
bar are probably contaminated. 



timation with normally distributed noise. For the arcs, a (mean) 
dust temperature of 87±7K (f3= 1.04±0. 12) is derived and for 
the bar 64±2 K(/3= 1.0±0.05) (see Fig.[l2£] For the arcs, 91% 
of the (3- values have values <1.1, and the maximum value for 
(3 is 1.88 with only 0.4% having values in the interval between 
1.75 and 1.90. For the bar, the maximum /3-value is 1.02. Low 
values for the spectral index have already been found in cir- 
cumstellar environments (e.g., |Knapp et al.||1993] l. Our result 
of a low spectral index (around 1) for quite high dust temper- 
atures around 65 K (bar) and 85 K (arcs) is in good agreement 
with the inverse temperature dependence of the dust spectral 



index as derived by Dupac et al. ( 2003) >. According to Dupac 



et al. 



( 2003| l, a low spectral index value might be explained by 
(1) the occurrence of large grains, (2) the fact that warm regions 
could harbor aggregates of silicates, porous graphite, or amor- 
phous carbon, having a spectral index around 1, or (3) an in- 
trinsic dependence of the spectral index on the temperature due 
to quantum processes, where in the case of high temperatures 



processes such as thermally-activated relaxation processes and 
temperature-dependent absorption could dominate. The chem- 
ical composition of the arcs might contain layered amorphous 
silicate grains (cfr. Verh oelst et al.|2006] > or larger grains (as re- 
cently detected around the AGB star W Hya, Norris et al.|2012| l 
explaining the low /3-index. In the case where the linear bar has 
an interstellar origin (see Sect. [5]), amorphous silicates — an im- 
portant component of the interstellar dust budget — might ex- 
plain the low /3-index as well. However, very small grains whose 
nature is thought to be carbonaceous ( |Desert et al.|19 86) might 
also be the origin of the low /3-index. Observations in the near-IR 
are prerequisites to constraining the detailed chemical content in 
the arcs and in the bar. 

A modified blackbody with a temperature around 63 K can 
explain the WISE 22 fim data of the linear bar, but not the WISE 
12 data. Admittedly, the flux in the WISE 12 fim band is 
quite uncertain and should be interpreted as an upper limit due 
to the large contribution of the central target (see Fig. [5]). It might 
be that the Herschel and AKARI images show the emission by 
larger dust grains compared to the WISE data, which can be fit- 
ted with a modified blackbody of ~ 145 K. Smaller grains indeed 
attain higher temperatures, hence emit at shorter wavelengths, 
than large grains, since the former absorb more efficiently per 
unit mass. 

The dust temperature map (for j3 = 1) and hydrogen column 
density^] in the arcs and linear bar as derived from the Herschel 



PACS data are shown in Fig. 13 A lower limit for the flux val- 
ues of 9, 8, and 4 MJy/sr was used for the 70, 100, and 160 fim 
data, respectively. Mainly in the eastern region, the flux values 
in the 160 image are too low to derive the dust temperature 
and hydrogen column density properly. More physical informa- 
tion on the eastern region can be obtained from either the slope 
of a linear fit to the flux values at 70,100, and 160 fim or from 
the ratio of the flux at 70 fim to the flux at 100/im (see bot- 
tom panels in Fig. 13 i.The arcs and bar have comparable dust 
colour temperatures. The dust temperature is higher both in the 
multiple arcs and in the linear bar in the direction of space mo- 
tion of Betelgeuse. A gradient in temperature is visible, with the 
temperature decreasing (and column density increasing) at larger 
distances from the central target. 

Considering energy conservation, one can show that the en- 
ergy flux of the gas entering the shock front per unit area and per 
unit time is around five orders of magnitude lower than the stel- 
lar radiation energy flux at the distance of the arcs and bar. As a 
result, the additional heating of the dust grain from the hot post- 
shock gas is negligible, and the temperature of the dust grain 
is dictated by radiative equilibrium conditions when only taking 
the stellar radiation into account (see also Sect. |4.2.T) . Assuming 
a smooth continuous outflow, the temperature of an iron or amor- 
phous silicate grain only heated by the stellar radiation is in the 
order of ranges from about 45 K to 60 K at a distance of 280- 
530" away from the central star (see Appendix [d]). The corre- 
spondence between the mean dust temperatures derived from the 
data and these theoretical calculations is satisfactory, taken the 
assumptions of the theoretical model into account. Specifically, 
one can imagine that the dust temperatures in a clumpy enve- 
lope might indeed be somewhat higher than the values shown in 
Fig. _ 



D.l 



since stellar photons can probably reach the outer en- 
velope regions more easily. In addition, smaller grains will attain 
a higher temperature since they absorb more efficiently per unit 
mass. 



The Herschel/SPIRE flux values are well off the fit, presumably due 
to contamination by ISM dust. 



3 for a dust opacity k„ = Kniis/vct) 13 , with uo = 1000 GHz, 
K = 0.1cm 2 /g, and/3=l JBeckwith et al.|l990 
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Figure 13: Dust temperature [K] (fi = 1, upper left), hydrogen column density [cm 2 ] (upper right), slope of the linear fit to the flux 
values at 70, 100, and 160 /im (bottom left), and ratio between the flux values at 70 and 100 (bottom right) as derived from the 
Herschel PACS images. See text for more details. 



3.2.4. Mass of the arcs and linear bar 

For an opacity, k\, of 185 c m 2 /g at 70 iim, 90 cm 2 /g at 100 ^m, 
and of 35 cm /g at 160 /im ( |Mennellaetal.|l998|l and assuming 
a dust-to-gas ratio of 0.002 ( Verhoelst et al.|2006) , we deduce 
a total gas and dust mass of ~2.4 x lO^M© for the arcs and 
^2.1 x 10 -3 Mq for the linear bar. The derived (gas+dust) mass 
is very sensitive to the assumed dust temperature and dust-to-gas 
ratio. Lowering the dust temperature in both the arcs and bar to 
30 K would yield a mass of 0.07 M Q and 0.029 M , respectively. 



Assuming that the linear bar has an ISM origin (see Sect|5]) with 
a dust-to-gas ratio around 0.01, the mass in the bar would be 
around 0.01 M . 

Assuming a smooth outflow, a first theoretical estimate for 
the mass in the bow shock shell can be derived from the equation 
dMohamed et al.|2012] > 



AL 



hell 



MR so 



(1) 
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with i?so the stand-off distance of the bow shock and v w the 
wind velocity. This expression assumes the mass is distributed 
over Att steradian in a spherical shell and only takes the contri- 
bution from the stellar wind into account. Thus for 9 ^ 90°, 
the expected shell mass would be half of Afghan • A bow shock 
shell mass of 2.4 x 10 -3 M Q for the arcs would imply a mass- 
loss rate of 2 x 1CP 7 M /yr (for a wind velocity of 14.5 km/s). 
However, CO and Hi observations tracing the inner envelope 
regions yield a mass-loss rate that is on e order of magnitude 
highe r ~2 x 10~ 6 M /yr (e. g., Paper II; |Rodgers & Glassgold 
|1991[ |Le Bertre et al.||2012| >. This hints that the mean mass- 
loss rate might have varied substantially, as also suggested in 
Sect. 



The hydrodynamical models presented in Sect. 4.2.1 yield a 
mass in the arcs of ~0. 1M e (including 90% C SM and 10% ISM 
material). According to Mohamed et al.]( |2012| l, the observed low 
mass in the arcs might indicate that the bow shock created by the 
RSG wind has a very young age, of the order of 20 000 yr, and 
may not yet have reached a steady state. But, the existence of 
(non-homogeneous) clumps in the inner wind travelling all the 
way to the bow shock region might also explain the low mass in 
the arcs as derived from the Herschel images. 

4. Discussion: Origin of multiple arcs 

In this section, we focus on the origin of the multiple arcs 
seen in the bow shock around Betelgeuse. Secti on|4.1| gives de- 
tails on some observational properties. In Sect. 4.2 the obser- 
vational data are compared with hydrodynamical simulations. In 
Sect. |4.3| we discuss the origin of the multiple arcs based on the 
observational and theoretical constraints. 



4. 1. Constraints from the observations 

From the observations described in previous sections, two con- 
cerns might be important for constraining the origin of the mul- 
tiple arcs seen in the bow shock region of Betelgeuse. (1) The 
different arcs are reasonably well fitted with concentric ellipses 
(see Fig. B.l in the Appendix |Bj). (2) The inner envelope shows 
clear evidence of a pronounced asymmetric clumpy structure. 
This leads to the hypothesis that the origin of the different arcs is 
a projection effect on the plane of the sky of the contact discon- 
tinuity created from the collision between different inhomoge- 
neous mass-loss events and the ISM. As shown in Fig. 6 of Cox 
et al. (2012 1, different inclination angles lead to projected out- 
lines of the contact discontinuity, where concentric ellipses with 
a smaller inclination angle yields a smaller projected size of the 
bow shock outline. This would mean that the arcs are shaped by 
an interaction of the inhomogeneous stellar mass loss at past and 
present epochs. 

Another explanation for the different arcs (and also linear 
bar) might be inferred from the recent results of |van M arie et al. 
( 201 l| l, who shows that small grains follow the movement of 
the gas mass elements, while larger CSM grains tend to pene- 
trate deeper into the ISM (see also Sect. |4.2[ ). This analysis takes 
only the inertia of the grains into account. However, since the 
grains are charged, they will gyrate around the magnetic field 
lines. Assuming a gas temperature in the shock of 10 000 K and 
a shock velocity equal to the space velocity, the typical Larmor 
radius for a grain with a size of 5 nm is ~ 1 x 10 14 cm, while 
larger grains with a size of 1 /im have a Larmor radius in the or- 
der of 4 x 10 18 cm. With a typical width of ~20" for the arcs, 
grains with a size JSO.I /im are position-coupled to the hot gas. 
One then could postulate that arcs further away from the central 



target could contain larger dust grains (^0. 1 /im), while the inner 
arcs contain the smallest grains coupled to the gas. In the case 
of randomized mass-loss variations (as suggested by the data 
tracing the inner envelope), one could consider the possibility 
that each mass-loss event has a favourable grain-size spectrum 
(i.e., the grain size distribution function reaches a maximum for 
some specific grain sizes). An analogous situation is obtained by 
Fleischer (1994), who modelled the dust formation in the case 
of a dynamical atmosphere and showed that far out in the wind, 
where the grain growth is definitely stopped, the grain-size spec- 
trum shows distinct grain-size peaks. A larger grain size implies 
a higher drift velocity and a larger Larmor radius, which eventu- 
ally might lead to a bow shock arc further away from the central 
target. However, it remains remarkable that only Betelgeuse in 
the whole MESS sample shows the appearance of multiple arcs 
in the bow shock region, while the same argument as given above 
might hold for other targets in the sample, especially the Mira- 
type variables with high-amplitude variations being favourable 
for a grain-size spectrum with strong distinct peaks. In addition, 
the dust temperature maps (Fig. 13 1 show no significant tem- 



perature difference between the different arcs, while it is well 
known that smaller dust grains are generally hotter than larger 
grains since the former absorb more efficiently per unit mass. 

The observed fragmentation of the bow shock might also be 
understood in terms of the effects predicted by Dgan i & Soker| 
(1998). They postulate that RT instabilities might fragment a 
bow shock in the direction of motion, enabling the ISM to flow 
into the inner parts of the envelope and creating a bow shock well 
inside the almost spherical, but very filamentary, haloes. For tar- 
gets close to the Galactic plane (as is the case of Betelgeuse, 
with b = —9°), the ISM magnetic field can make RT instabil- 
ties very efficient. This might explain the recent detection by 
|Le Bertre et al.| ( |2012| l of a detached H I shell with a radius of 
2'. Magnetic fields would suppress certain modes but accentu- 
ate others, changing the appearance of the envelope, and poten- 
tially leading to 'RT rolls' (or stripes). An example resembling 
the image of Betelgeuse in shown in Fig. 3a in Dgani & Soker 
( 1998). We should remark, however, that the calculations as per- 
formed by Dgani & Soker ( 1998 ) were for the case of a planetary 
nebula, which have much higher stellar wind velocities than red 
supergiants or AGB stars. In the sample of 78 AGB stars and 
supergiants analysed by |Cox et al. ( 2012| ), Betelgeuse is the tar- 
get closest to the Galactic plane (z = — 5pc), and it is the only 
target showing this multiple arc-like structure. Four other targets 
in the MESS sample are at similar distances from the Galactic 
plane (z < 20 pc), one being a supergiant (/i Cep) which re- 
sembles Betelgeuse in different aspects, including a similar dust 



mass-loss rate to Betelgeuse (Josselin & Plez 2 007| ). The main 
difference for parameters influencing the bow shock morphol- 
ogy is the wind velo city, which is a factor 2.4 higher in fi Cep 
( peBeck et al.|2010| . While for Betel geuse v±/v w > 1, this ra- 
tio is < 1 for n Cep in which case it has been predicted by Dgani 
et al.|([T996]) that the bow shock is stabler. The Herschel images 



presented by |Cox et al. ( 2012[ ) show, that (possible) instabilities 
in the bow shock of /x Cep have comparable sizes to Betelgeuse, 
while different arcs are indeed not seen. 

If present, the strength of the ISM magnetic field can be pre- 
dicted from the angular separation between the arcs. An angular 
separation of 30" at a distance of 300" yields an Alfven speed 
in the pre-shocked ISM of ^4 km/s (Eq. 4 in|Dgani|1 998). For 
an ISM density of 1.9 cm " 3 ( |Cox et al.|20i2] ) and assuming that 



mainly H I is present in the bow shock, this yields a magnetic 
field of 3 /zG. The average strength of the total magnetic field 
in the Milky Way is about 6 /iG near the Sun (Biermann 2001) 
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and increases to 20-40 fiG in the Galactic centre region. Dense 
clouds of cold molecular gas can host fields of up to several mG 
strength ( |Heiles & Crutcher|200"5] ). 

4.2. Constraints from hydrodynamical modelling 

For insight into the origin of the complex CSM-ISM interface 
and the non-occurrence of large-scale instabilities in the arcs, 
we have simulated the hydrodynamical evolution of the enve- 
lope of Betelgeuse. The simulations also show if dust grains are 
good tracers for the gas-driven dynamics in the bow shock. The 
hydrodynamical computations were done using MPI-AMRVAC 
code (Keppens et al. 2012| van Marie et al. 201 1) . More de- 
tails on the numerical method and initial conditions are given in 
Appendix [E] 

Table 2: Physical parameters of the envelope of Betelgeuse. 



Parameter 



Value 



Photospheric temperature 
Stellar radius 
Stellar mass 
Mass loss rate 
Dust-to-gas ratio 
Distance 

Terminal wind velociy 
Space velocity w.r.t. ISM 



T* = 3500K (a) 
i?*=7.5xl0 13 cm (a) 
M* = 20M Q (i,) 
M = 2.5 x lCT 6 M Q /yr (c) 
V> = 0.002 (d) 
D=197 pc (e) 
Woo = 14.5 km/s (c) 
i;* = 28.3km/s (/) 





Rodgers & Glassgold 


« 


1991); (b) Dolan et al. 


(2008); (c) pa 


per II; 


(d) 


Verhoelst et al. (20061; 


e) Harper et al. 1 2008}; {f > Ueta et al. 


20081 



The simulation was run for five different situations, using 
variations of a basic model to explore several aspects of the bow 
shock morphology and its dependence on the properties of the 
wind and the ambient medium. The physical parameters for the 
basic model are given in Table|2] The gas density in the ambient 
medium is set at 3 x 10~ 24 g cm~ 3 (or rin ~ 2 cm~ 3 ) and a 
temperature of 100 K (reflecting the temperature of a cold neu- 
tral medium). The five different simulations have the following 
specifications: 

- The first simulation (A) uses the basic wind parameters that 
are specified in Table [2] and the ambient medium specified 
above. We assume that all dust grains have the same size, 
5 nm, and density, 3.3 g cm~ 3 . 

- The second simulation (B) is the same as simulation A, ex- 
cept for the size of the dust grains, which is increased to 
100 nm to demonstrate the behaviour of large dust grains. 

- The third simulation (C) has a warm ambient medium of 
8 000 K, reflecting the temperature of a warm neutral or par- 
tially ionized medium. The ambient medium is assumed to 
be heated by an outside radiation source. Therefore, whereas 
the ambient medium is kept at a minimum temperature 
of 8 000 K, the stellar wind, which is protected from UV- 
radiation by the ionized gas of the bow shock, has a mini- 
mum temperature of 100 K. 

- The fourth simulation (D) has the same ambient medium 
temperature as simulation C, but with lower ambient medium 
density (5 x 10 -25 gem -3 ) and higher space velocity (u* = 
72.5 km/s). This simulation was added to compare our re- 
sults with the recent hydrodynamical models of Mo hamed] 
|et~aL| ( |20T2l >. 

- For the fifth simulation (E), as for the first and second, we set 
the minimum temperature to 100 K throughout the entire do- 



main. However, instead of a smooth wind, we assume a peri- 
odic variation in the mass loss rate. For 1 000 years the mass- 
loss rate is assumed to be high (M = 2.5 x 10~ 5 M Q /yr), then, 
for the next 9000 years, it is low (M = 2.5 x 10~ 8 M /yr). 
This 'picket-fence' type variation is then repeated period- 
ically. The two mass-loss rates are related as A/high = 
1000 Af low The total amount of mass lost over the 10000 
year period is the same as for the first two simulations. The 
period of 10 000 years reflects the typical timescale for AGB 
thermal pulses. 

The effects of changing one or more input parameters in the 
simulations are discussed in Sects. 14.2. 11 - 14.2.31 A link to the 
online movies is provided in the online Appendix[F] A reflection 
on the origin of the multiple arcs based on the outcome of the 
simulations is given in Sect. |4.3| 

4.2.1 . Cold ambient medium 

The result for simulation A with an ambient medium temper- 
ature of 100K is shown in Fig. [14] After a simulation time of 
300 000 yr, the place of the bow shock interaction has stabi- 
lized: the termination shocl^] occurs at 0.25 pc (or 262") and 
the bowshoclj^] at 0.34 pc (or ^356"), in good agreement with 
the observations (see Sect. 3.2 neglecting the inclination an- 
gle). The turbulent astropause/astrosheatfr] thus has a width of 
^90". Although both the wind termination shock and the for- 
ward shock are smooth, the contact discontinuity, which sep- 
arates the shocked wind from the shocked ambient medium, 
shows instabilities. These do not penetrate the shocks on either 
side, but cause deviations in the shape of the forward shock. 



Figure 14 shows three types of instabilities in the gas. (1) A 
thin layer of shocked wind material, right on the contact discon- 
tinuity, has a higher density and lower temperature than the sur- 
rounding gas. This is the result of the radiative cooling instabil- 
ity, which causes dense gas to cool faster. As a result, dense gas 
loses thermal pressure and is compressed, which increases the 
density and therefore the cooling rate. (2) The high density gas 
of the shocked wind extends "fingers" into the shocked ambient 
medium, which form in the forward shock and increase in size 
as they travel downstream. These are Rayleigh-Taylor instabili- 
ties, resulting from low density, but high-pressure, material (the 
shocked ambient medium) pushing back denser shocked wind. 
(3) Due to the shear-force between the shocked wind and the 
shocked ambient medium, the Rayleigh-Taylor instabilities de- 
form and show a circular motion. These are Kelvin-Helmholtz 
instabilities and are primarily visible downstream from the star. 
All three features were also found by van Marie et al. ( 2011| l. 

The temperatures in the shocked gas are lower than the 
kinetic energy of the free-streaming wind and the ambient 
medium. For a fully adiabatic shock, the temperatures at the 
front of the bow shock would be about 27 300 K for the shocked 
interstellar medium and around 7 200 K for the shocked wind. 
Instead, the maximum temperature for the shocked interstellar 
medium is 14 000 K and the shocked wind does not reach more 
than 4 000 K. [van Marie et al. ( |2011| l find similar temperatures. 
This is mainly due to the effective radiative cooling. Also, as 
the thermalized gas moves downstream, it decompresses, trading 



4 location where the wind velocity goes from supersonic to subsonic 
values 

5 location where the ISM velocity goes from supersonic to subsonic 
values 

6 location where the ISM pressure equals the CSM pressure 
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Figure 14: The dust particle density [cm 3 ] and gas density [gem -3 ] (left); and the gas velocity [km/s] and gas temperature [K] 
(right) for a simulation with an ambient medium temperature of 100 K (simulation A). Although the shocks are smooth, the contact 
discontinuity shows extensive instabilities, which start small at the front of the shock and grow as they move downstream. The dust 
concentration follows the contact discontinuity. 



temperature for an increasing velocity, which causes the tail of 
the shock to be cooler than the front. This can be seen by com- 
paring the absolute velocity and temperature of the gas (right 
panel of Fig. 14 1. 

Comparing the gas and dust density (left panel of Fig. \\A\ 
shows that the dust particles penetrate the shocked wind layer, 
but are brought to a stop at the contact discontinuity, which they 
tend to follow. In that sense, small dust grains are good tracers of 
the gas-driven dynamics at the contact discontinuity. This con- 
firms the earlier results found by van Marie et al. (2011 1. As a 
result, the instabilities at the contact discontinuity should also be 
visible in the infrared (assuming that the observations have suffi- 
cient resolution, which is the case for the Herschel instruments). 
Since the dust is concentrated at the contact discontinuity, other 
features that are typical of the shocks would most likely not show 
up in infrared images. For example, in the region behind the star, 
the wind termination shock curves back toward the axis of mo- 
tion as its thermal pressure balances against the ram pressure of 
the wind. This effect is completely invisible in the dust. The be- 
haviour of the circumstellar environment over the full 200000 
year period, starting from the moment when radiative cooling is 
introduced, is presented online in animated form (Appendix |F}. 
This movie shows that the number of instabilities increases with 
the introduction of radiative cooling. These instabilities form in 
the front of the bow shock and then travel downstream. 

The total mass of the gas trapped inside the bow shock re- 
gion, found by integrating the density between the termination 
shock and bow shock for < 9 < 90°, is approximately 0.1 M Q 
of which 90% is shocked wind material. The shell mass is con- 
siderably higher than the value obtained from the Herschel data 



(see Sect. 3.2 1. As already suggested, a possible explanation may 



be that the envelope around Betelgeuse is very clumpy or that 
Betelgeuse is a very young RSG star so that the bow shock is 
not yet fully formed (Moham ed et al.|2012 i. This could also ex- 
plain the smoothness of the bow shock, since instabilities need 



time to form. Most of the mass is situated at the sides of the 
star, where the volume of the shock region is much larger; e.g., 
when integrated only over < < 45° the total mass becomes 
1.8 x 10~ 2 M s , of which about one third is shocked ISM. 

Assuming radiative-equilibrium balance, the heating of the 
dust grain due to energy transfer from the hot post-shock gas 
can be estimated. An approximation is given by Eq. (5.40) from 
Tielens (2005), which assumes the absorption efficiency of the 
dust as a function of wavelength to be a simple power law: 
Q(A) = Qo(A/Ao)' 3 , with Q Q = 1, A = 2na and /3 = 1. As a 
result, the dust temperature only depends on the mean intensity 
of the radiation field J and the size of the dust grains a, 



Td oc 



0.2 



(2) 



The total energy flux entering the shock front is 



total 



= F, 



wind 



2 Airr 2 



ISM 



1 



:PISM 



5.5 x 10 5 erg cm 



(3) 



Assuming that only 20% of the energy flux is converted to radi- 
ation, the additional heating from the dust grain due to the hot 
post-shock gas is only ^6 K, while the heating due to the stellar 
radiation is one order of magnitude higher (see Fig. |D.l| in the 
Appendix |Dj). 

The effect of grain size on the distribution of dust in the bow 
shock region is demonstrated in Fig. [15] This shows the result 
for simulation B for a larger grain with size 100 nm. Instead 
of following the contact discontinuity, the larger dust grains are 
concentrated at the forward shock and, in the tail of the bow 
shock, even penetrates the undisturbed ambient medium. This 
behaviour, which was also shown in |van Marie e t al. (20lT|), is 
the result of the larger momentum of large grains, compared to 
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Figure 15: Similar to Fig 14 but for a simulation with large (100 nm) dust grains (simulation B). The dust distribution is completely 
different, tracing the forward shock, rather than the contact discontinuity, and even penetrating the unshocked interstellar medium 
(see text for more details). 



the drag force. The large grains are less tightly coupled to the gas 
and can therefore penetrate deeper into the bow shock region 
and the ambient medium beyond. However, since larger grains 
are less numerous than smaller ones, the infrared morphology 
of the bow shock will mainly be determined by the smaller dust 
species. 



4.2.2. Warm ambient medium 

When the ambient medium is kept at a minimum temperature 
of 8 000 K (simulation C), the results look different, as shown 
in Fig. 



The result of simulation D with the same high ambient 
medium temperature, but lower ambient medium density and 



16 The Rayleigh-Taylor instabilities are much smaller. 



This is the result of radiative cooling. The instabilities, which are 
made up of stellar wind material, are allowed to cool to 100K, 
whereas the surrounding shocked ambient medium has a mini- 
mum temperature of 8 000 K. Therefore, the gas in the Rayleigh- 
Taylor fingers has a lower thermal pressure than the surrounding 
gas and is compressed, inhibiting their growth. This is different 
from the results shown in |Cox et al.| ( |2012| ), which also showed 
large instabilities for the warm ambient medium. However, these 
may have been partially due to numerical effects of modelling a 
spherical wind expansion on a cylindrical grid, which lead to a 
large instability close to the symmetry axis. 

The layer of shocked interstellar matter extends further from 
the axis of motion and has an almost triangular shape. The com- 
bination of the smaller instabilities and the larger distance be- 
tween contact discontinuity and forward shock also reduces the 
effect of the instabilities on the shape of the shock, which is com- 
pletely unaffected. 

Because of the high thermal pressure of the ambient medium, 
the wind region is more constrained; in fact, the free-streaming 
wind has not yet reached the lower boundary, but is still con- 
strained by the shocked wind. This was also seen in the results 
of |Cox et al.| l |2012) . As in the first simulation, the small dust 
grains follow the contact discontinuity. 



higher space velocity, is shown in Fig. 17 The high stellar veloc- 
ity creates a stronger shear force between the shocked wind and 
the shocked ambient medium. As a result, the Kelvin-Helmholtz 
instabilities are much greater, leading to a large-scale instability 
that changes the shape of the entire bow shock. This is contrary 
to the results of Moham ed et al.] (|2012 ), who obtained a smooth 
bow shock for similar input parameters. This difference may be 
the result of a different numerical treatment (grid versus SPH), 
but may also result from the total timescale of the simulation, 
because Kelvin-Helmholtz instabilities need time to grow (see 
accompanying movie in Appendix |F|. Owing to the larger col- 
lision speed at the bow shock, the shock temperatures are much 
higher (~ 10 5 K) than for the other simulations. This largely 
negates the effect of the high ambient medium temperature, since 
its thermal energy is much lower than the kinetic energy. While 
the general structure of the bow shock is unstable, the local dust 
shell directly in front of the star is quite smooth. Also, the in- 
stability of the bow shock leads to a structure that curves in the 
opposite direction from the wind termination shock. Similar fea- 
tures have been found by Wareing et al. (2007a) and Cox et al. 
| |2012) . 

4.2.3. Clumpy mass-loss episodes 

The hypothesis that the different arcs might be the result of the 
collision of different clumpy mass-loss episodes having filling 
factors <1 with the surrounding ISM (see Sect. 



4. 1 1 can be 



checked using the hydrodynamical simulations. In the first in- 
stance, this can be done in a simplified 2D scheme, where the 
collision between the ISM and a CSE with a variable mass-loss 
rate has been modelled. A full 3D model is beyond the scope 
of this paper. As outlined above, we simulated this effect by 
varying the mass-loss rate by a factor 1000 over a period of 
10 000 yr. Since the change in dust density at the edge of a clump 
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Figure 16: Similar to Fig 14 but with an ambient medium temperature of 8 000 K (simulation C). The instabilities are much smaller. 
Since the shocked interstellar matter does not cool below 8 000 K, the shocked ambient medium layer extends further outward than 
in Fig. 14 On the other hand, the shocked wind region is more restricted. 
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Figure 17: Similar to Figs. 14 and 15 but with an ambient medium temperature of 8 000K, an ambient medium density of 
10 -25 gem -3 , and a stellar velocity of 72.5 km/s (simulation D). The bow shock is both locally and globally unstable, the shock 
temperature of the bow shock is high, because of the higher collision speed. 

is very steep, the mass-loss rate variations were implemented in the shells hit the wind termination shock before they can expand 
a 'picket-fence' way. far enough. 

The variation in mass-loss rate (and therefore in ram pres- 
sure) causes a widespread disturbance in the morphology of 
the bow shock region. The wind termination shock, where the 
isotropic thermal pressure of the shocked gas is equal to the ram 
pressure of the wind, shows a wave-like pattern in the downwind 
region as it conforms to the variable ram pressure of the wind. 
At the front of the bow shock the variation in the free-streaming 



Figure 18 shows the result of simulation E with variable mass 



loss. The mass leaves the star concentrated in shells that are sep- 
arated by voids with a density that is three orders of magnitude 
lower. As they travel outward, the shells expand owing to their 
own internal pressure. In the wake of the star the shells eventu- 
ally merge. To the front of the star this does not occur because 
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Figure 19: Zoom of Fig. 18 The white arrow indicates the possibility of multiple arcs in the front of the bow shock region. 



wind increases the instability at the contact discontinuity. The 
resulting Rayleigh-Taylor and Kelvin-Helmholtz instabilities are 
more numerous than in the case of a smooth wind (Fig. 14 1 and 
show more small-scale structure. 

Finally, the impact of subsequent shell collisions means 
that multiple arc-like structures appear between the termina- 
tion shock and the contact discontinuity (for a detailed view see 
Fig.[T9|. This may be an indication of the origin of the multiple 
arcs observed in the bow shock of Betelgeuse. 

The animation of the simulation (Appendix|F|i shows the in- 
fluence of the variable mass-loss rate, in terms of both the for- 
mation and growth of local features, as well as of the large-scale 



morphology of the bow shock. Despite the variation in ram pres- 
sure, the location of the bow shock mostly remains unchanged. 
This is the result of the period of the variation relative to the dy- 
namical timescale of the bow shock region. The bow shock loca- 
tion can only change if the dynamical timescale is significantly 
shorter than the period of the variation. This is not the case in 
this simulation. The bow shock region has a cross-section of the 
bow shock of approximately 0. 1 pc at the front and up to a par- 
sec to the sides of the star. With a sound speed of about lOkm/s 
(for a temperature of 10 000 K), the timescale varies between 
10000 and 100000 years. Since the period of the variations is 



18 



L. Decin et al.: The enigmatic envelope of Betelgeuse 



only about 10000 years, the bow shock does not have enough 
time to adjust to the variations. 

4.3. Discussion 

A. Multiple arcs: The hydrodynamical simulations show some 
good agreement between the observations and the theoretical 
predictions, which might point us toward the origin of the arcs; 
e.g., the observed ratio between the projected — minimal — dis- 
tance between the star and the bow shock outline in the direction 
of relative motion and the distance perpendicular to that, i.e. at 
9 = 90°, is —0.66 for the different arcs. The same ratio is ob- 
tained in simulations A, C, and E, while the ratio is significantly 
lower for simulations B and D. 

The hydrodynamical simulations show that the formation of 
dust clumps in the inner wind regions, simulated in all its gen- 
erality by implying spherically symmetric mass-loss variations, 
can have an impact on the observed morphology of the bow 
shock region, and potentially might lead to the formation of mul- 
tiple arcs in the bow shock. Dust clumps with higher grain size 
have a larger drift velocity, large Larmor radius, and longer stop- 
ping length when passing through the bow shock, and eventually 
might create a bow shock arc further away from the central tar- 
get. However, as said above, we would then expect that other tar- 
gets in the MESS sample would show multiple arcs in the bow 
shock region as well. 

An effect not included in the MPI-AMRVAC-code is photo- 
ionization. In the case of the CO waves (or ripples) detected on 
the surface of the Orion molecular cloud, Ber ne et al.| (2010) 
argues that ultraviolet radiation has created a (small) insulating 
photo-ablative layer, allowing the development of a KH insta- 
bility with a wavelength at least one order of magnitude longer 
than the insulating layer. A speculative idea is that the multiple 
arcs found around Betelgeuse are undulations created by photo- 
chemical effects stimulating the growth of some instability with 
wavelength around the observed width of the arcs, which are 
then traced by the smaller dust grains that tend to follow the gas 
instabilities. However, also for this effect, it remains surprising 
that Betelgeuse is the only target in the MESS sample showing 
multiple arcs. 

In our opinion, the most probably reason for the origin of 
the multiple arcs is the combined effect of (1) a turbulent non- 
homogeneous mass loss, which might enable the ISM to flow 
into the envelope and create multiple bow shocks, together with 
(2) the effect of the interstellar magnetic field, which can make 
some instabilities very efficient while suppressing other modes 
(see Sect. 4.1 1. The angular separation between the arcs points 
toward a magnetic field of —3 fj,G. 



B. Suppression of large-scale instabilities: Remarkably, the 
size of (local) density clumps tracing potential instabilities in- 
side the arcs as seen in the Herschel images (wi th an upper limit 
to the maximum length of —30", see Sect. 3.2 1 is often smaller 



-30", see Sect. 

than predicted by the hydrodynamical simulations. Most of the 
instabilities seen in Figs. 14 and 19 are Rayleigh-Taylor (RT) in- 
stabilities, which are slightly deformed by the action of Kelvin- 
Helmholtz (KH) instabilities. Thin-shell instabilities such as the 
NTSI ( |Vishniac|1994l l or TAI ( |Dgani et al.|1996| l do not occur, 
since the shell between both shocks is too thick: with a typical 
thickness of the shell in the order of 0. 1 pc, the crossing time is 
in the order of 10 4 yr (for a velocity of — lOkm/s) by which time 
an individual mass element is already displaced by a few times 
10 17 cm, preventing the growth of an instability. Thermal insta- 



bilities are also ruled out, since radiative shocks with a power- 
law cooling functi on (A tx T 7 ) are stable against thermal insta- 
bilities if 7 > 1 ( |Gaetzet al.|1988 1. 

The only example in our simulations showing only small- 
scale instabilities is the case of a warm ambient medium around 
Betelgeuse (with a temperature of —8000 K) for which a strong 
growth of the instabilities is prevented by compressing the 
Rayleigh-Taylor fingers. As shown in Fig. 16 instabilities in 



the direction of the space motion of the target then have typi- 
cal lengths < 1 x 10 17 cm (or 34"). 

We should also realize that the size and growth rate of 
the instabilities are directly related to the density of the gas 
( Chandrasekhar 1961}. In our simulations, compression due to 
radiative cooling increases the density of the shocked wind. 
However, it is possible that we over-estimate the cooling. The en- 
ergy loss due to radiative cooling is proportional to the ion den- 
sity times the electron density. If the gas is only partially ionized 
(as in the 3 000-8 000 K regime), the cooling may be reduced, 
which in turn would reduce the compression and therefore the 
density contrast at the contact discontinuity, damping the insta- 
bilities. 

Not considered in our simulations is the effect of a mag- 
netic field. Depending on the direction of the magnetic field with 
respect to the shear flow, a magnetic field might (de)stabilize 
KH instabilities (|Miura & Pritchett|T982l |Keppens et al.|1999) . 
Miur a & Pritchett| ( |1982[ ) found that in a sheared magnetohydro- 
dynamical flow in a compressible plasma, modes with fcA < 2 
are unstable, with A the scale length of the shear layer and k the 
wavenumber. Compressibility and a magnetic field component 
parallel to the flow (Bo || vn) are found to be stabilizing effects. 
For the transverse case (Bo -1 Vo), only the fast magnetosonic 
mode is destabilized, but if k ■ Bo 7^ 0, the instability contains 
Alfven-mode and slow-mode components as well. Keppen s~et al.| 
(1999) studied the case of an initial magnetic field aligned with 
the shear flow. In the case the initial magnetic field is unidirec- 
tional everywhere (uniform case), the growth of KH instabilities 
is compressed, while a reversed field (i.e. a field that changes 
sign at the interphase) acts to destabilize the linear phase of the 
KH instability. As a result, the non-occurrence of large-scale in- 
stabilities might potentially be explained by the presence of a 
magnetic field. As suggested in Sect. 4. 1 the detection of multi- 
ple arcs in the bow shock might also point towards the presence 
of an ISM magnetic field. With a strength of —3 ^,G, the mag- 
netic pressure is ~ 3.6 x 10~ 13 dyn cm~ 2 , being two orders of 
magnitude lower than the ram pressure at the contact disconti- 
nuity. Without dedicated simulations, it is not possible to predict 
how the presence of a weak ISM magnetic field might influence 
the growth rate of the instabilities. 

One might also argue that the bow shock shell is still very 
young so that instabilities might not have had enough time to 
grow. Using the prescriptions of Brighenti & D'Ercole (1995) 
the typical growth time for RT instabilities to reach a size of 30 
at a distance of —200 pc is only —7000 yr, while the KH growth 
timescale is at least an order of magnitude smalle^] However, 
according to our simulations, the shell would not have had time 
yet to reach the observed stand-off distance. Changing the mass- 
loss rate by a factor of — 10 only reduces the current stabilization 
time from - 1 00 000 yr to -30 000 yr. 

The clumpy structures found in the inner envelope region 
might also have an imprint on the appearance of the bow shock. 



However, one should realize that for KH instabilities the difference 
in velocity Av > 2« soun( j, with Av being at maximum — 43km/s for 
these simulations. 
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With shock velocities below 50km/s, the percentage of dust 
mass destroyed (i.e. returned to the gas phase) is negligible 
(Tielens 2005 ). But, one should also take into account that dust 
clumps tend to dissipate. One possible depletion mechanism is 
heat conduction, which is thought to be unimportant in the case 
of cool supergiant (and AGB) winds owing to the small dif- 
ference in temperature between the clump and its environment. 
According to |Hartquist et al. ( 1986| l, clouds are, however, ab- 
lated when they are in relative motion with respect to the sur- 
rounding medium. Following |Smith et al.| ( fl988] l, we estimate 
that all dust clumps with a size smaller than ~ 10 15 cm (or 0.3" 
at a distance of 200 pc) will be ablated by the time the clump 
reaches the bow shock region; larger dust clumps might still be 
visible in the bow shock region. In the case where most of the 
dust grains formed in the inner wind had already been ablated by 
the time they arrived at the bow shock region, the clumps seen in 
the outer arcs might also reflect ISM inhomogeneities progres- 
sively engulfed by the expanding shell. 

5. Discussion: Origin of the linear bar 

It is striking that the arcs and the linear bar show the same far- 
infrared colour temperature and density (as deduced from the 
PACS data). This might suggest that they are illuminated by the 
same source, probably Betelgeuse itself, and belong to the same 
physical phenomenon (i.e. the bow shock). As already argued in 
previous section, larger grains are capable of leaving the shocked 
wind layer, and can penetrate the shocked or even unshocked 
ISM. The large grains almost completely ignore the morphology 
of the shocked gas shell and penetrate the unshocked ISM a con- 
siderable distance (~0. 1-0.2 pc, [van Marie et al. 2011). Since 
they are very low in number, they will not influence the gas mor- 
phology. However, even for very large grains, the shape of the 
bow shock still shows some curvature (see Fig. 4 in |van M arie 
et al.|20lT), w hich is not seen in the Herschel images. Recently, 
Mackey et al. (2012) have suggested that the linear bar is the relic 



of a blue supergiant (BSG) wind interaction with the ISM, where 
Betelgeuse has recently entered the red supergiant (RSG) phase. 
The BSG reverse shock collapses in on the RSG wind producing 
an inner shocked shell (the arcs in the Herschel image). The r~ 2 
decrease in dust emission makes the contact discontinuity look 
a little more like a bar (compared to the projected density in a 
2D simulation) because the strongly curved parts at larger radii 
are fainter. If so, Mackey et aL]( |2012[ l predict that the bar should 
eventually begin to curve back into a bow shock shape, with a 
radius of curvature that can be larger than the current stand-off 
distance of the bar suggests. The large FOV Herschel images 
cannot give extra support to this proposed scenario since no cur- 
vature is seen. 

We thus should consider the possibility that the linear bar is 
not related to Betelgeuse, but is an interstellar structure, which 
is by accident co-spaced with Betelgeuse. If so, the motion of 
the star and angular separation between the star and bar imply a 
collision between the two in about 17500 yr, while the arcs and 
linear bar would collide in about 5000 yr. 

The linear bar might be the edge of an interstellar cloud. The 
results obtained by Dic key & Garwood] ([1989 ) prove that it is 
not unreasonable that Betelgeuse is near an interstellar cloud. 
Dickey & Garwood ( 1989 ) showed that a few hundred clouds are 



present within 1 kpc having a column density around 10 cm 



an extended circumstellar wake associated with the CSE around 
the AGB star X Her. Mapping shows that the star lies (in pro- 
jection) near the periphery of a much larger H I cloud that also 
exhibits signatures of interaction with the ISM. While the radial 
velocity of the cloud near X Her overlaps that of the CSE, this is 
not the case for a Ori (see Figs.[7]l. However, the probability of 
a mere chance superposition in position, velocity, and direction 
of space motion between a cloud and target is extremely low. 
Betelgeuse might thus be close to an interstellar cloud, but have 
a different radial velocity. 

The bar might also be a linear filament in the interstellar cir- 
rus, as also seen in the Galactic centre, whose possible origin is 
linked to the Galactic magnetic field ( |Yusef-Zadeh et al~[ [l984; 
|Yusef-Zadeh & Morris|1987 ). Linear filaments might be caused 
by particle acceleration by shocks related to preexisting filamen- 
tary magnetic structure in the ISM ( |Heyvaerts et al.| 19 88 ) or by 
particle energization at reconnection sites at the ionized surfaces 
of molecular clouds and subsequent particle loading of ISM 
magnetic field lines (Serabyn & Morris 1994). Filamentation can 
also be the result of synchrotron thermal instabilites ( |Bodo et IT 
1990): both the thermal condensation mode and slow magne- 
tosonic wave can create structures that are aligned with the back- 



ground magnetic field. |Rosner & Bodo ( 1996 1 show that a par- 



ticle acceleration process, akin to the acceleration of the anoma- 
lous cosmic -ray component associated with the solar wind ter- 
mination shock, assisted by radiative instabilities driven by syn- 
chrotron emission can provide an explanation for the filamenta- 
tion process. Some fraction of accelerated particles will be line- 
tied to the external (ISM) magnetic field. The local field mag- 
netic field is parallel to the Galactic plane and follows the local 
spiral arms dAndreasyan & Makar ov|1989| ), having a pitch angle 
of about —8° ( |Han|2008[ ). Optical polarization measurements in 
the direction of Betelgeuse on a scale of several degrees show 
that the position angle of the local ISM galactic field (|Bin gham 
|& Shakeshaftl[T967l |Heiles||2000| |Heiles & Crutcher||2005| ) is 
consistent with the position of the linear bar. However, typical 
temperatures of interstellar cirrus clouds are around 15 K (e.g. 

2010| ), much lower than the tempera- 



Miville-Deschenes et al. 



(cfr. Fig. [13) . With an estimated mass for the linear bar around 
^0.01 M , the molecula r cloud mass spectrum has a density of 
~10pc~ 3 (see Fig . 2 in|Dickey & Garwood||l989[ ). Using Hi 
21 cm observations, Mat thews et al.| ( |201 \\ has recently detected 



tures we deduced from the Herschel maps. Only when the fila- 
mentary cirrus is by chance close to Betelgeuse can the radiation 
of the target illuminate and heat the filamentary structure. 

6. Conclusions 

This paper presented new Herschel PACS and SPIRE images of 
the red supergiant Betelgeuse. The Herschel data show that the 
collision between the wind of the red supergiant and the inter- 
stellar medium has created a very complex structure. Multiple 
arcs at a distance of 6-7' are detected with typical dust temper- 
atures in the range of 40-140 K. At a distance of ~9' from the 
central target, an intriguing linear bar is detected, with the same 
colour temperature as the arcs. Betelgeuse is the only target in 
the MESS sample showing these two phenomena. 

The inner envelope structure is very inhomogeneous and tes- 
tifies to a turbulent mass-loss history. The Herschel images show 
the first evidence of a high degree of non-homogeneous distri- 
bution for the material beyond 15", which even persists until the 
material collides with the ISM. This inhomogeneity is proba- 
bly linked to the giant convection cells recently detected in the 
outer atmosphere of Betelgeuse, which might enable localized 
dust creation and ejection. 

The Herschel data have been complemented with data from 
the UV (from GALEX), near-IR (from WISE), and radio (from 
the GALFA-HI survey) to disentangle the different components 
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that might result in this complex bow-shock structure and in the 
linear bar. The average temperature in the arcs is ^85 K and in 
the bar is ~63 K. A low dust emissivity index, j3, around 1 is 
deduced from fitting the spectral energy distribution. This might 
point toward (1) the occurrence of large grains, (2) a grain com- 
position resembling layered amorphous silicates grains, porous 
graphite, or carbonaceous grains, or (3) high temperature pro- 
cesses influencing the quantum effects. The dust temperature 
in the arcs and bar is regulated by the stellar radiation field of 
Betelgeuse, while heating due to the hot post-shock gas is negli- 
gible. 

The Herschel PACS images of the inner envelope and the H I 
results by |Le Bertre et al.| ( |2012| l suggest that the mean dust and 
gas density might have changed considerably ~32 000yr ago 
(radius of ~2'). The idea of a strong variation in the mean den- 
sity is also supported by the low bow shock shell mass deduced 
from the Herschel data. This change in mean density might be 
from a change in mean mass-loss rate or the formation of an in- 
ner bow shock well inside a filamentary halo, in which case the 
ISM is able to flow into the inner parts of the envelope creating 
a bow shock. 

Detailed hydrodynamical simulations were performed to un- 
derstand the morphology of the bow shock around Betelgeuse. 
The stand-off distance and the ratio between the project distance 
between the star and the bow shock outline in the direction of rel- 
ative motion and the distance perpendicular to that are predicted 
well by the hydrodynamical simulations, except in the case of 
large dust grains (around 100 nm; simulation B) or when a low 
ambient medium density and high space velocity are considered 
(simulation D). It is clear that small dust grains are good tracers 
for the gas-driven dynamics at the contact discontinuity. Most 
of the mass in the bow shock is situated to the sides of the star, 
and only ^20% is situated in a cone with angular extent <45°. 
The dust temperature in the bow shock is mainly determined by 
heating by the stellar radiation, while energy transfer from the 
hot post-shock gas can be neglected. Increasing the temperature 
in either the ambient interstellar medium or in the shocked wind 
prevents the growth of strong instabilities. Allowing for strong 
variations in the mass-loss rate in the inner envelope has an im- 
pact on the appearance of the bow shock, and eventually multiple 
arc -like structures might appear between the termination shock 
and the contact discontinuity. 

Based on the observations and on hydrodynamical simula- 
tions, different hypotheses are formulated to explain the origin 
of the multiple arcs and to understand why no large-scale in- 
stabilities are seen in the bow shock region. In our opinion, the 
two main ingredients to explain both features are (1) a clumpy 
mass-loss process and (2) the influence of the Galactic magnetic 
field. The occurrence of clumps might also explain the low to- 
tal (gas and dust) mass in the arcs as deduced from the data 
(~ 2.5 x 10~ 3 M Q ). The Galactic magnetic field might serve 
as an extra trigger to fragment the bow shock region. The an- 
gular separation between the arcs seen in the Herschel data is 
compatible with a magnetic field of ^3 fiG, in agreement with 
the average strength of the magnetic field in the Milky Way. 

The linear bar might be the edge of an interstellar cloud il- 
luminated by Betelgeuse or a linear filament whose a possible 
origin is linked to the Galactic magnetic field. Since no curva- 
ture is present in the bar, we believe that the bar is not directly 
linked to a previous blue supergiant wind, as recently discussed 
by |Mackey etaT1 ( (20T2) >. 
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Appendix A: Previous images of the bow shock 
around Betelgeuse 

For completeness, we show in this section the bow shock around 



Betelgeuse as observed with IRAS by Noriega-Crespo et al. 
( fT997l > (Fig. [AH} and AKARI by |Ueta et al.| ( |2008| l (Fig.|A3|~ 
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Figure A. 2: AKARI/FIS false-colour maps of a Ori in the SW 
bands — N60 (65 /zm) and WIDE-S (90 //m) at 15'7pixel scale 
— and in the LW bands — WIDE-L (140 /im) and N160 
(160/Ltm) at 30"/pixel scale — from left to right, respectively 
(Ueta et al. 2008) >. Background emission has been subtracted 
by a combination of temporal filters during data reduction. RA 
and Dec offsets (with respect to the stellar peak) are given in ar- 
cminutes. The wedges at the top indicate the log scale of surface 
brightness in MJy sr -1 . North is up, and east to the left. 



Figure A. 1: IRAS 60 /jm image of the bow shock around 
Betelgeuse as detected bylNoriega-Crespo et al. ( 1997 1. 



Appendix B: Additional figures 

Fig. |B.1| PACS 70 /im image on which concentric ellipses were 
fitted to the different arcs. 




Appendix C: Herschel PACS image of R Leo 

The Herschel PACS 70 /mi image of the AGB star R Leo are 
compared to a hydrodynamical simulation in Fig. C. 1 Rayleigh- 
Taylor instabilities, slightly deformed owing to the action of 
Kelvin-Helmholtz instabilities, are visible both in the data and 
in the simulations. 



Appendix D: Dust grain temperature in the 
envelope around Betelgeuse 

Under the assumption of a smooth continuous outflow and us- 
ing the parameters as specified in Table |2j the dust grain tem- 
peratures for Fe and amorphous silicates in the envelope around 
Betelgeuse have been calculated using the MCMAX code (Min 
|et al.|[2009] ). Only the stellar radiation is taken into account as 
the energy source to compute the dust temperatures assuming 
radiative equilibrium. The particle shape is represented by CDE 
(continuous distribution of ellipsoids Bohren & Huffman 1983) 
for particles of which the equivalent spherical grain would have 
a grain size of 0.1 /xm. The composition of the amorphous sili- 
cates is taken from |de Vries e t al. (2010), who derived that the 




Figure B.l: PACS 70 /im image overplotted in red with three 
concentric ellipses with the centre at a = 5h54m56.56s and 
<5 = 7°22'16.18" radii of 622", 585", and 546" (as measured 
from the centre) and for a position angle of 47.7°. 



infrared spectrum of the AGB star Mira is best fitted with 65.7% 
MgFeSiC>4 and 34.3% Mg 2 SiC»4, implying an average composi- 
tion of Mgi.36Feo.64SiC<4, having a /3-index of ~1.8. The result- 
ing dust temperatures are shown in Fig. |D.1| The temperature 
for the amorphous silicate grains is lower than for the Fe-grains, 
owing to the higher opacity of Fe at the near-IR wavelengths, 
where the star emits most of its photons. 
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Figur ed: PACS 70 fjm image of R Leo (left |Cox etaT 
|2012| l compared to a hydrodynamical simulation of the wind- 
ISM interaction computed with the AMRVAC code (van Marie 



et al. 



201 l| l. The parameters for this simulation are a stel- 



lar wind velocity of 15km/s, a constant gas mass-loss rate of 
1 x 10~ 6 M Q /yr, a dust-to-gas mass ratio of 0.01, a space veloc- 
ity of 25km/s, a local ISM density of 2 cm' 3 , and an ISM tem- 
perature of 10K. The upper right figure shows the gas density, 
which ranges (on log scale) between 10 -24 and 10 _19 g/cm 3 ; 
the lower right figure represents the dust grain particle density, 



ranging (on log scale) between 10 10 and 10 3,0 cm 
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Figure D.l: The dust-temperature profiles in the envelope of 
Betelgeuse as modeled with MCMAX ( |Min et alT]|2009| ). The 
different lines indicate the dust temperature of a specific dust 
species: red for amorphous silicates and blue for metallic iron. 
The full black line gives the mean dust temperature profile, as- 
suming thermal contact between iron and amorphous silicate 
grains. The green lines show a power-law distribution for the 
dust temperature for j3=l. The vertical dashed line indicates the 
inner radius of the dust shell. 



Appendix E: Numerical simulations of the bow 
shock 

E. 1. Numerical method 

For the simulations the MP I -AMRVAC hydrodynamics code 
(Keppens et al. 2012|l was used. F or a comparison to similar 
work, we refer to | van Marie et al.| ( [201 1) . The MP I -AMRVAC 
code solves the conservation equations for mass, momentum, 
and energy on an adaptive mesh grid. It includes optically 
thin radiative cooling, using the exact integration algorithm 



( |Townsend]|2009[ ), and a cooling curve for gas at solar metal- 
licity (Wang Ye, private communication). For these simulations, 
a two-dimensional approach was taken, assuming an inclination 
angle of 0°. As deduced by Cox et al. (2012 1, the inclination an- 
gle of the space motion with respect to the plane of the sky is 
small (-8°). 

Dust was included in the simulations using a two-fluid 
approach, where the dust is treated as a pressure-less fluid 
( Paardeko oper & Mellema|2006[ |Woitke|2006a| |van Marie et al.| 
|201 l[|Cox et al.|2012|i. The dust is linked to the gas through the 
drag force ( Kwok|1975] ): 



fd = (1 - oi T ) 7r n d p a AvJ(Av) 



(E.l) 



with rid the dust particle density, p the gas density, a the radius of 
the dust grain, vt = y/ikT/nih the thermal speed, and Av the 
velocity difference between dust and gas. For a given mass, the 
particle density decreases with a~ 3 . Therefore, the drag force is 
weaker for larger dust grains. The sticking coefficient for gas- 
dust collisions, ax, is estimated as 



a T = 0.1 + 0.35 e 



V( T s+ T d)/ 500 



(E.2) 



(Dec in et al.|2006[ l, with Td the dust temperature and T g the gas 
temperature. This is an improvement over our earlier simulations 
( |van Marie et al.|201 1[ Cox et al. 2012[ l, where we kept the stick- 
ing coefficient constant. Because the dust temperature is typi- 
cally much lower than the gas temperature in the shocked CSM- 
ISM region, we approximate the sticking coefficient by setting 
T d = C0 

Unlike our previous models (van M arie et al.|201 l[|Cox et al.| 
|2012| ), which used a cylindrical grid, these simulations are run 
on a two-dimensional spherical grid with a physical domain of 
R = 3pc by 9 = 180°. The grid has a minimal resolution of 
160 radial by 80 azimuthal grid cells and five additional levels 
of refinement are allowed, giving the simulations an effective 
resolution of 2 560 x 1 280 grid cells. 

E.2. Initial conditions 



At the inner radial boundary, the wind parameters are set ac- 
cording to the values in Table [2] Since the star is moving rel- 
ative to the ambient medium, the outer radial boundary is di- 
vided into two zones: (1) for 9 < 90° an inflow boundary is set 
with the ambient gas entering the grid as a parallel stream with 
v = 28.3 km/s, (2) for 9 > 90° the gas is allowed to flow out o f 
the grid. This is similar to the approach of |Villaver et al. (2012). 
The gas density in the ambient medium is set at 3 x lO - ^ g 
cm~ 3 (or riH ~ 2 cm" 3 ) and a temperature of 100 K (reflect- 
ing the temperature of a cold neutral medium). This temperature 
also serves as a lower limit for the gas temperature throughout 
the physical domain. We assume that there is no dust in the am- 
bient medium. 



This approximation is not valid in the region inside the termination 
shock where gas and dust temperatures might be similar far away from 
the star. As described by |van Marie et al.|p0lT) , this region is not sim- 
ulated in our computations: the wind interaction is initialized by filling 
a spherical area of radius 0.1 pc around the origin with wind material 
having equal gas and dust velocity. As seen in Fig. [O] the observed 
dust temperature in the wind-ISM interaction region is ~70 K, which 
would yield an increase of the sticking coefficient of maximum 2.8% 
compared to the results of Eq. |E.2| for Td = 0. In the case of large dust 
grains which penetrate a cold ISM (simulation B), the change in the 
sticking coefficient in this ISM region would be ~10%. 
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To avoid numerical instabilities we start the simulation with- 
out radiative cooling and let it run for 100 000 years. By that 
time the bow shock has reached its equilibrium position (where 
the ram pressure of the wind is balanced by the ram pressure of 
motion of the star through the ambient medium). From this point 
in time we include radiative cooling and let the simulation run 
for 200 000 years. 
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Appendix F: Animations of bow shock 
simulations 

Figure F. 1 : Movie of AorLTISMlOOK showing gas density and dust particle density for the 200 000 year period from the introduc- 
tion of radiative cooling to the end of the simulation for the model with a cold ambient medium. 



Figure F.2: Similar to Fig. F. 1 AorLTISM8000K showing the 200 000 year period from the introduction of radiative cooling to the 
end of the simulation for the model with a warm ambient medium. 



Figure F.3: Similar to Fig. Fl the third movie (AorLvardM) shows the 200000 year period from the introduction of radiative 
cooling to the end of the simulation for the model with a cold ambient medium and a varying mass-loss rate. 



Figure F.4: Similar to Fig. F.2 the fourth movie (AorLTISM8000K_highvel) shows the 200000 year period from the introduction 
of radiative cooling to the end of the simulation for the model with a warm ambient medium and a high stellar velocity. 



Figure F.5: Similar to Fig. F. 1 the fifth movie (AoriJargegrain) shows the 200000 year period from the introduction of radiative 



cooling to the end of the simulation for the model with a cool ambient medium, but for large (100 nm) dust grains. 
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